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C/3 . ABSTRACT 

Cold dark matter models predict the existence of a large number of substructures 
within dark matter halos. If the cold dark matter consists of weakly interacting mas- 
sive particles, their annihilation within these substructures could lead to diffuse GeV 
' emission that would dominate over the annihilation signal of the host halo. In this work 

Q\ , we search for GeV emission from three nearby galaxy clusters: Coma, Virgo and For- 

nax. We first remove known extragalactic and galactic diffuse gamma-ray backgrounds 
r*"** " and point sources from the Fermi 2-year catalog and find a significant residual diffuse 

emission in all three clusters. We then investigate whether this emission is due to 
■ (i) unresolved point sources; (ii) dark matter annihilation; or (iii) cosmic rays (CR). 

f^) Using 45 months of Fermi-LAT data we detect several new point sources (not present 

fSJ . in the Fermi 2-year point source catalogue) which contaminate the signal previously 

analyzed by Han et al. Including these and accounting for the effects of undetected 
point sources, we find no significant detection of extended emission from the three 
clusters studied. Instead, we determine upper limits on emission due to dark matter 
' annihilation and cosmic rays. For Fornax and Virgo the limits on CR emission are 

consistent with theoretical models, but for Coma the upper limit is a factor of 2 below 
the theoretical expectation. Allowing for systematic uncertainties associated with the 
treatment of CR, the upper limits on the cross section for dark matter annihilation 
from our clusters are more stringent than those from analyses of dwarf galaxies in the 
Milky Way. We rule out the thermal cross section for supersymmetric dark matter 
particles for masses as large as 100 GeV (depending on the annihilation channel). 
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1 INTRODUCTION 

The existence of dark matter (DM) in the universe has so 
far only been deduced from its gravitational effect, due to 
the lack of electromagnetic interactions of the DM with 
itself or with baryonic matter. There are several elemen- 
tary particle candidates for DM in various extensions o f 
the standard model of particle physics l|Bertone et al.|[2004 ). 
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Weakly interacting massive particles or WIMPs (with a self- 
interaction cross-section at roughly the weak scale) are one 
of the popular dark matter candidates. These particles could 
be related to the electroweak symmetry breaking which 
is currently being explored by experiments at the LHC. 
For example, within the framework of the minimal super- 
symmetric standard model (MSSM), the lightest neutralino 
emerges as a candidate WIMP that is stable over cosmologi- 
cal timescales and can annihilate into standard model parti- 
cles. WIMPs behave as cold dark matter since their primor- 
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dial velocity dispersion is negligible. High resolution N-body 
simulations show that cold dark matter halos contain a pop- 
ulation of self-bound substructures (subhalos) whose num- 
ber decreases w ith increasing subhal o mass as N oc M~ a 
with a » 1.9 jDiemand et all 120071 ; ISpringel et ail |200S| ; 
iGao et al.ll201lh 

Much effort has been devoted to the search for WIMPs 
either directly or indirectly. Direct detection involves iden- 
tifying the rare events of DM scattering off ordinary mat- 
ter or searching for new particles near the weak scale at 
the LHC. Indirect detection involves looking for the anni- 
hilation or decay products of dark matter in cosmic rays 
and gamma rays. In particular, pair annihilation produces 
gamma-ray photons at a rate proportional to the square 
of the dark matter density, which then propagate, almost 
without absorption, to the observer. In this case, the Galac- 
tic c enter should be the brightest gamma-ray source on the 
sky ijSpringel et al.ll2008l . and references therein). Extended 
emission (distinct from the central point source) was re- 
ported from the centra l 1° ar ound the Galactic center by 
iHooper fc Goodenougbl l|201lf ): see also IHooper fe Linden! 
|2011bl ). This emission has bee n interpreted as a signal from 
dark matter annihilation, but iBovarskv et al.l |201 if ) point 
out that, since its radial profile is consistent that of known 
Fermi point sources, the emission could all originate from 
a point source at the Galactic centerQ An intriguing as- 
pect of a DM explanation for the gamma-ray emission from 
the Galactic centre is that the inferred particle mass of 
around 10 GeV is also the mass claimed to be required to 
explain other data, such as the s ynchrotron emission from 
the Milky Way's ra dio filament (jLinden et al. l201lh and 
the "WMAP Haze" [Finkbeinerl | 2004| ; IHooper et all [20071 ; 
IHooper fc Lindenl l2011al ). as well as signa ls from the di- 
rect detection exp eriments DAMA/LIB RA (|Bernabei et alj 
l2010h. CoGeNT dAalseth et all l2011blH) and CRESST-II 
( Angloher et~aD l201lT l. These signals, however, are incon- 
sistent with other direct d etection experime nts, such as 
CDM S \ Ahmed et al.l l201lh and Xenon- 10 (|Aprile et all 
l201ll ). We refer the reader to iHooperl l|2012h for review. 

It has recently been reported that the 7-ray emission 
from the region around the Galacti c center exhibits a line- 
like excess at energies ~ 130 GeV dBringmann et a"1l2012l : 
IWenigedlioij iTempel et afll2012i ; ISu fe Finkbeinerl l2012h . 
The intepretation of this signal as arising from dark mat- 
ter p articles, however, is controversial (see IBovarskv et alJ 
120121 ). 

Targeting the entire sky rather than the Galactic cen- 
ter in searching for annihilation radiation may seem a good 
strategy since this takes advantage of the large-scale dis- 
tribution of dark matter while avoiding some of the uncer- 
tainties arising from the astrophysical modelling of galac- 
tic gamma-ray sources. However, the fact that we are lo- 
cated near the center of the Galactic halo and that most of 
the annihilation emission outside the Galactic center is pro- 
duced by dark matter substructures (|Diemand et al.l [20071 ; 



1 See also a prelim inary result by the Fermi-LAT collaboration 
llVitale et al.|[2009h ■ which is consistent with either DM annihila- 
tion from a halo with a density profile, cuspier than the typical 
NFW col d dark matter galactic ha lo (an inner slope of ~ —1.25 
to -1.4) l|Navarro et al.ll 19961. Il997h . or with proton collisions ac- 
celerated by the central galactic black hole. 



ISpringel et al.| [2008) results in a gamma-ray map from anni- 
hilation that is almost uniform on large scales. This makes 
detection within the Milky Way halo a difficult task, exacer- 
bated by the additional uncertainty of having to model the 
extrag alactic background, which is more important on large 
scales l|Zahariias et al.ll20ld : lBaxter fc Dodelsonll201ll ). 

Dwarf galaxies are the most DM-dominated ob- 
jects known, are relatively free from astrophysical con- 
tamination and appear compact on the sky. They are 
therefore promising targets to search for DM anni- 
hilation radiation. Recent joint analyses of eight to 
ten dwarf galaxies ( Gcringcr-Samcth & Koushiappas l201ll : 
iThe Fermi-LAT Collaboration: M. Ackermann et al.l 120111 ) 



resulted in no significant detection but have began to 
rule out the canonical annihilation cross-section of 3 x 
10 _26 cm 3 s _1 for DM masses below ~ 30 - 40 GeV. 

Galaxy clusters are the most massive virialized DM 
structures in the universe and are also good targets for indi- 
rect DM searches. The presence of a large population of DM 
substructures (or subhalos) predicted by numerical simula- 
tions further enhances the detectability of DM in clusters. 
Although the total mass within subhalos amounts to only 
10 to 20 percent of the total halo mass, the density enhance- 
ment within subhalos can boost the total cluster annihila- 
tion luminosity by a factor as high as 1000 when extrapo- 
lated down to a subhalo mass limit of one Earth mass, the 
fiducial cutoff in the primoridal power spe ctrum of density 
fluctuations for a typical 100 GeV WIMP (|Gao et al.ll201ll ; 
iPinzke et aill201ll 'l. As the distribution of subhalos is much 
less concentrated than that of the smooth main halo, the 
total annihilation emission from clusters is predicted to be 
extended. Thus, attempts to detect DM annihilation assum- 
ing a point source or NFW-squared profile could miss most 
of the signal. In fact, just such a search using the 11-month 
Fermi-LAT data has yi elded no significnat dete ction of emis- 
sion from six clusters ll Ackermann et aT1l2010r i. 

Using the 45-month data, we consider possible contri- 
butions from cosmic ray (CR) induced gamma-ray emission 
and from DM annihilation. For the former (which can be 
as high as, or higher than the emission from cluster DM 
annihilation dJeltema et al 1 l2009l ; IPinzke fc PfrommeijlioTol ; 
IPinzke et aTH201ir). we adopt the semi- analytic method de- 
veloped by IPinzke fc Pfrommerl feOlCh. For the later, we 
adopt the model proposed by IGao et all (|20111 ) for the 
cluster DM annihilation profile. We provide constraints on 
both the CR and DM components for the three galaxy 
clusters analyzed by Han et al. (2011): Coma, which is 
predicted to hav e the highest signal-to-noise according to 
IGao et al.ll|201ll l, and Fornax and Virgo which are predicted 
to have the lowest a strophysical contamination according to 
IPinzke et ail (|201lh . 

The current paper replaces an earlier version by a subset 
of the authors (arXiv:1201.1003). After submission of that 
version, it was pointed out to us that a number of point 
sources are present in the full three-year LAT data which 
were not detected significantly in the data used for the "of- 
ficial" Fermi point source catalogue available at the time of 
our analysis, the LAT 2-year p o int so urce catalogue (2FGL; 
IThe Fermi-LAT Collaboration! (|201lh 'l . We now carry out 
our own point source detection in the regions of interest 
and find several new point sources. We account for these 
new detections in our analysis, as well as for the fact that 
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a significant part of the "smooth" extragalactic background 
is contributed by point sources below the detection thresh- 
old; this alters the noise properties of this background. Both 
changes reduce the significance of the diffuse components ap- 
parently detected in the first version of our paper, so that 
we can now reliably onl y place upper limits. 

iHuang et al.l ()201 lh have recently reported a failure to 
detect significant DM annihilation emission from a com- 
bined analysis of eight galaxy clusters. Our work differs from 
theirs in several respects: firstly, we assume a DM annihi- 
lation profile based on high resolution cosmological simula- 
tions (|Gao et al.l l201ll ); secondly, we assess the impact of 
cosmic rays in the detection of dark matter; and finally, we 
include in our sample the Virgo cluster which turns out to 
be the best candidate. The constraints we set o n the annihi- 
lation cross-section are consistent with those of lHuang et al.l 
|201l| y 

The paper is organized as follows. In section [T] we de- 
scribe the data and provide an overview of the models of 
the Virgo, Fornax and Coma galaxy clusters regions used 
in the analysis (see Table [T]). The specification of the non- 
standard components of the models (dark matter and cosmic 
rays brightness profiles) is provided in Sec. [2] The constrains 
on CR emission and DM annihilation that we obtain are 
summarized in section [3] and discussed in Sec. [4] 

The cosmological para meters used in thi s work are the 
same as those assumed by iGao et al.l (|201ll ): fi m = 0.25, 
Qa = 0.75, h = 0.73. 



1.1 Data preparation 

We analyze the first 45 months of data (04/08/2008 to 
20/05/2012) from the Fermi-L AT, trimmed with the cuts 
listed below, to select high quality photon events. This typ- 
ically results in ~ 10 5 photons within a radius of 10 de- 
grees around each cluster, while the expected number of 
annihilation photons is of the order of 10 2 according to 
Fig. [3] The most recent instrument response function (IRF) , 
P7CLEAN_V6, is adopted for the analysis, in accordance 
with our event class selection^ The resulting gamma-ray 
images for the three clusters are shown in the left panel of 
Fig. Q] for Virgo and in Fig. IC II for Coma and Fornax. 



Minimum Energy 
Maximum Energy 
Maximum zenith angleS 
Event Class0 

DATA-QUAlFI 

LAT CONFIGJ 

ABS (ROCK ANGLE jf| 

ROI-based zenith angle cut 



100 MeV 
100 GeV 
100 degrees 
3 (P7CLEAN) 
1 
1 

< 52 degrees 
yes 



2 http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi 

3 We also tried using P7SOURCE.V6 IRF and Event Class 2 
data. The results are consistent with those presented in this pa- 
per. 



We list the basic properties of the three clusters in Ta- 
ble [1] 



1.2 Maximum-likelihood fitting 

We use the pyLikelihood tool shipped with the Fermi 
Science Tools software package (version v9r27pl-fssc- 
20120410) to perfor m a maximum likelihood (ML) analysis 
l|Mattox et al.l ll996). After applying appropriate data cuts, 
as described in section [LTI we bin the data into 0.1 degree- 
wide pixels and 30 logarithmic energy bins within a radius 
of 10 degrees around each cluster. This large radius is cho- 
sen to account for the large LAT PSF size at low energies 
(4 ~ 10 degrees at 100 MeVQ). An exposure cube is com- 
puted around each cluster covering 25 degrees in radius and 
the 30 energy bins, using the gtexpcube2 tool. 

In the standard Fermi likelihood analysis, the photon 
counts within each pixel are treated assuming Poisson statis- 
tics for each energy bin to calculate the likelihood. The best- 
fit parameters are obtained when the likelihood for the entire 
data set is maximized. The significance of a given component 
of interest (e.g. DM or CR) from the ML fitting is quantified 
by the likelihood ratio statistic, 



TS = -2]n(L /L), 



(1) 



where L is the maximum likelihood for the full model and 
Lo is the maximum likelihood for the null hypothesis, i.e, 
the model without the component of interest. According to 
Wilk's theorem, this test statistic, TS, approximately fol- 
lows a x 2 distribution when the null hypothesis is true, with 
one degree of freedom for our case where the normalization 
is the only extra parameter in the alternative model. The 
probability that a given value of TS arises purely from fluc- 
tuations of the null hypothesis is: 



■dx. 



(2) 



The factor — comes from the constraint that the normaliza- 
tion parameter be non-negative. The significance of a detec- 
tion can thus be quoted as y/TSa (one-sided Gaussian con- 
fidence). Upper limits on the extra normalization parameter 
N are obtained by searching for a null hypothesis L' where 
N in the full model is constrained to be equal to the upper 



5 ZENITH ANGLE (degrees): angle between the reconstructed 
event direction and the zenith line (originates at the center of the 
Earth and passes through the center of mass of the spacecraft, 
pointing outward). The Earth's limb lies at a zenith angle of 113 
degrees. 

6 EVENT CLASS: flag indicating the probability of the event 
being a photon and the quality of the event reconstruction. 

7 DATA-QUAL: flag indicating the quality of the LAT data, 
where 1 = OK, 2 = waiting review, 3 = good with bad parts, 
= bad 

8 LAT-CONFIG: flag for the configuration of the lat (1 = nomi- 
nal science configuration, = not recommended for analysis) 

9 ROCK ANGLE: angle of the spacecraft z-axis from the zenith 
(positive values indicate a rock toward the north). 

9 The LAT PSF size scales roughly as E~°- a , so at 1 GeV it is 
~ ldeg 
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Table 1. Basic Properties of Target Clusters 





CoiTlel 


J_ Ul llClA. 




RA (deg) 


194.9468 


54.6686 


187.6958 


DEC (deg) 


27.9388 


-35.3103 


12.3369 


D A (Mpc) a 


95.8 


17.5 


16.8 


M 200 (M ) b 


1.3el5 


2.4el4 


7.5el4 


f200 (deg) b 


1.3 


4.1 


6.2 


Jnfw c 


5.9e-5 


4.1e-4 


1.2e-3 


Enhancement due to subhalos within r200 


1.3e3 


6.5e2 


1.0e3 



"Angular diameter distance, from the NASA extragalactic database for Coma and Fornax, and from lTullv &: Sh ava (1983) f° r Virgo. 
'Cluster halo mass defined as the mass within the radius, ^ 200 , within which th e average density equals 200 times the c ritical density o f 
the universe. Values for Coma and Fornax are taken from [Pinzke et al. 

1 11201111 . while the value for Virgo is taken from lTullv fc Shaval 

||1984| 1. 

integrated coefficient, Ji n t = f^n Jdfl, over the solid angle spanned by the cluster virial radius, assuming a smooth NFW density 

profile. 

^Enhancement to the total annihilation luminosity within the virial radius due to substructures, extrapolated to a subhalo mass limit of 

io- 6 m q 



limit, Nul, so that ln(L /L) — —1.35, corresponding to the 
95% confidence interval. 

1.3 Model 

For the analysis we constructed a model to fit the data in- 
cluding all known foreground and background emission, as 
well as DM and CR components, as appropriate. We in- 
clude all the point sources from 2FGL within a radius of 
15 degrees from the cluster center in the model, plus the 
most recent galactic (GAL) and extragalactic (EG) diffuse 
emission given by the template files gal_2yearp7v6_v0.f its 
and iso_p7v6clean.txt. Additionally, we have searched the 
45-month data for new point sources; we detect several of 
them within the cluster region (see Appendix [A] for more 
detail) and these are also included in our model. The nor- 
malization of the GAL and EG diffuse components are al- 
lowed to vary during the fitting. Within the cluster virial 
radius there are two 2FGL point sources and one newly de- 
tected point sou rce in Fornax, six 2FGL, including the cen- 
tral AGN (M87: lAbdo et alJlioolh . plus four newly detected 
ones in Virgo. We allow the normalization and power-law 
spectral index of these thirteen point sources to vary freely. 
In addition, the parameters of all sources with variability in- 
dex greater than 50 located within 10 degrees of the cluster 
centres are allowed to vary. Parameters for the other point 
sources are fixed as in the 2FGL catalog. From now on we re- 
fer to the model with GAL, EG and the known point sources 
as the "base model". 

A DM annihilation surface brightness template (given 
by the dimensionless factor J, see Eqn.|4]in Sec. 12. 1[) is gen- 
erated for each cluster out to a 15 degree radius by summing 
up both the contribution from a smooth NFW profile and 
the contribution from subhalos. This J-map is used to fit for 
extended cluster annihilation emission. For the point source 
model, the integrated factor Jnfw (see Eqn. [5| is used to 
derive an annihilation cross-section from the fitted total flux. 



Similarly, a CR photon template is generated for each clus- 
ter out to three times the cluster virial radius, where the 
surface brightness has dropped to below 10~ 5 of the central 
value and beyond which the model is not reliable. Images for 
various model components are shown in Fig. [T]takmg Virgo 
as an example. We discuss these templates in more details 
in Sec. H 

In the traditional Fermi analysis, the EG template is 
treated as a smooth component where all emission below 
the nominal point source detection limit is assumed to come 
from a smoothly distributed diffuse component. In this work, 
we also consider a more realistic one where a fraction is 
assumed to be contributed by fainter point sources with a 
number-flux relation which extrapolates smoothly from that 
measured for brighter sources. In this case the photon counts 
within a given pixel are no longer Poisson-distributed since 
the photons arrive in packets. In principle, one can use the 
full distribution of photon counts from a population of ran- 
domly placed point sources to calculate the likelihoods L 
and Lo, but Eqns. [T] and [2] and the corresponding discus- 
sion, are not affected. However , sinc e the full distribution 
of photon counts in this case ijHanl I in prep,j l) is compli- 
cated and difficult to implement in the likelihood analysis, 
instead of recalculating L and Lo, in this work we use Monte- 
Carlo simulations to re-evaluate the distribution of TS for 
the more realistic background model and provide corrections 
to the results of the standard analysis where needed. 



2 MODELING GAMMA-RAY EMISSION IN 
CLUSTERS 

We model the observed gamma-ray emission in clusters with 
several components as shown in Fig. [T] the galactic fore- 
ground (GAL), the extragalactic background (EG), emis- 
sion from known point sources, DM annihilation and CR- 
induced emission. The GAL and EG diffuse emission are 
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Virgo 



DEC (deg) 




Figure 1. Decomposition of the Fermi-LAT image in the region of the Virgo cluster into model components. The observed photon count 
image from 100 McV to 100 GeV is shown on the left. The right panels show the integrated image over the same energy range for the 
various model components: galactic diffuse emission, extragalactic diffuse emission, detected point sources, cosmic-ray photons and DM 
annihilation emission, as labeled. The green dashed circle in the "Point Sources" panel marks the virial radius of the cluster. The small 
circles mark the newly detected point sources which are not present in the 2FGL. The "Fluctuation" panel shows the residual image for 
our best-fit DM model. The images have been enhanced individually in colour space for contrast. Note the apparent structure in the 
extragalactic component which is due to different exposure times at different positions. 



given by the most recent templates, gal_2yearp7v6_v0.f its 
and iso_p7v6clean.txt, which can be obtained from the 
Fermi-LAT data server, while the point sources include 
those from the LAT 2-year po i nt so urce catalogue, 2FGL 
ijThe Fermi-LAT Collaboration! l201ll ). as well as several, 
newly detected by us, in the 45-month data. In addition, 
an improved EG model which includes a population of un- 
detected sources is also analyzed. We now describe in detail 
our models for DM annihilation and CR emission. 



2.1 Dark matter annihilation emission 

The gamma-ray intensity along the line-of-sight due to DM 
annihilation is given by: 



8 7T 



dE 



< a fV > 



ij^fm, (3) 



where M x is the DM particle mass; p x the density of DM; 

jf the particle model dependent term giving the differ- 
dE 

ential number of photons produced from each annihilation 
event as a function of energy, E, in a particular annihila- 
tion channel, /; and < a/v > is the velocity-averaged cross- 
section (or annihilation rate) for that channel, which is pre- 
dicted to be constant in the low velocity li mit appropriate 
to pre sent-day cold DM particles (see e.g.. lJungman et ah! 
l|l996l )). The line-of-sight integration of the density squared 
is often expressed in terms of a dimensionless factor, 



J = 



:(; 



p x (l)dl. 



(4) 



8.5kpc v 0.3GeV/cm 3 Jlos 

If the source size is much smaller than the instrumental 
beam size, a point source approximation is applicable. In 



this case, the integration of J over a large enough solid an- 
gle, Af2, is used to determine the total flux for the point 
source, J int = / AfJ JdQ. 

The cluster annihilation e mission is modele d with the 
extended profile suggested by iGao et al.l (|201lh . However, 
for part of the analysis and for comparison purposes, we 
will also use the point source approximation which, although 
inappropriate, has been employed in all previous analysis of 
Fermi-LAT data from clusters. We shall refer to models that 
assume these two profiles respectively as EXT and PT. If the 
cluster follows a smooth NFW profile, then its integrated J 
factor which determines the total annihilation flux can be 
found as 

4i 2 '! 1 1 / 1 , 2 ,., 

3 D A 8.5kpc O.SGeV/cm 3 

Here Da is the angular diameter distance to the clus- 
ter and p„ and r 3 axe the characteristic density and ra- 
dius of the NFW profile. They are related to halo con- 
centration, c, and virial radius through the relations, p s — 

200 c 3 p c , . 

— r jj— — r- and r a = 7-200 /c, with p c the cnt- 

3 ln(l + c) — c/(l + c) 

ical density of the universe, r2oo the cluster virial radius 
within which the average density is 200p c and the con- 
centration parameter, c, is given by the following mass- 
concentration relation: 

M20O \ -0.097 



5.74( 



2 x lO^/i-iM, 







(6) 



(Duff v"et al.ll2008h . Here, M200 is the mass enclosed within 
r2oo- Extrapolating to a cutoff mass of 10 _6 Mq, the exis- 
tence of subhalos will increase this flux by a factor 



b(M 200 ) = Jsub/JNFW = 1.6 x lO" 3 (M 2O o/M ) 



(7) 



IGao et all (|201lh . Using the results of the simulations by 



6 J. Han et al. 



these authors, the surface brightness profile of subhalo emis- 
sion can be fitted within r2oo by the following formula: 



J S ub{r) 



m(M 20 o)jNFW 



7rln(17) 



+ 16r 2 



(r^raoo). (8) 



Below we fit the subhalo emission surface brightness beyond 
the virial radius and extrapolate to several times the virial 
radius using an exponential decay, 

Jsub(r) = Jsub(r20o)e- 2 - 37nr/r20 °- 1) (r > r 2 oo). (9) 

The total annihilation profile is the sum of the contributions 
from a smooth NFW profile and the subhalo emission. This 
is completely dominated by subhalo emission except in the 
very center of the cluster. We show the total annihilation 
profile and its decomposition into main halo and subhalo 
contributions in the left panel of Fig. [3] taking Virgo as an 
example. This profile is further inflated after convolution 
with the LAT point spread function. 

We consider three representative annihilation channels, 
namely into b—b, — and r + — r~ final states. The anni- 
hilation spectrum is c alculated using the DarkSUSY package 
jGondolo fc SiUJl999h , which tabulates simulation results 
from PYTHIA0 We also include the contribution from in- 
verse Compton (IC) scattered photons by energetic electron- 
positron pairs produced during the annihilation p r ocess , fol- 
lowing the procedure described in IPinzke et al.l |201lf ) . In 
general, three external energy sources are involved in the 
dissipation and scattering of the injected electrons from an- 
nihilation: the Cosmic Microwave Background (CMB), in- 
frared to UV light from stars and dus t, and the interstella r 
magnetic field. However, as shown bv lPinzke et al.l |201lf ). 
the latter two components are expected to be important 
only in the inner region of clusters (< 0.03r2oo), correspond- 
ing to less than 0.2 degrees for our three clusters. Including 
them would introduce a position-dependent component to 
the annihilation spectrum, so for simplicity we only consider 
the contribution of CMB photons in the IC calculation. For 
the bb channel, IC photons only contribute significantly to 
the low energy spectrum for relatively high neutralino mass, 
while for the leptonic channels, which have plenty of ener- 
getic electrons, the IC emission can completely dominate the 
annihilation emission over the full energy range of interest 
for the highest neutralino masses considered. 

We not e that the elec t rowea k correcti ons recently 
propo sed by ICiafaloni et al.l l|201ll ) (see also ICirelli et al.l 
(|20 1 lh ") can bring visible differences to the leptonic channel 
spectra at high WIMP masses before IC scattering. However, 
since IC photons dominate at the high mass end and the elec- 
troweak correction only significantly changes the positron 
yields at low energy, thus having little effect on the IC spec- 
trum, the electroweak correction to the total spectrum is 
still negligible. The total photon yields are shown in Fig. [2] 
The almost flat spectrum with a cutoff around the energy 
corresponding to the WIMP mass comes from prompt an- 
nihilation emission including continuum secondary photons 
and final state radiation from charged final state particles. 
The low energy rise originates from IC scattered CMB pho- 
tons. 



2.2 Cosmic-ray induced gamma-ray emission 
within clusters 

The cosmic ray induced gamma-ray emission is calculated 
following a semi-analytic prescription, derived from high res- 
olution numerical simulations of galax y clusters, that mod- 
els co smic ray physics self consistently l|Pinzke fc Pfrommerl 
2010). The gamma-ray photon production rate (or source 
function) from pion decay is found to be separable into a 
spatial and a spectral part: 



qcR(r, E) 



dN-, 



dtdVdE 



= A(r)s(E), 



(10) 



where the spatial part, A(r), is proportional to the square of 
the gas density profile multiplied by a slowly varying radial 
function parametrized by cluster mass. The spectral part, 
s(E), is almost independent of cluster mass and has a power- 
law form, dA r 7 /dln(_B 7 ) oc E~ 13 , for the energy range 1 ~ 
100 GeV but flattens at low energies, as shown in Fig. [2] 
We summarize the detailed form of A(r) and s(E) plus the 
gas density profile for the three clusters derived from X-ray 
observations in the Appendix. 

The differential gamma-ray flux from this source func- 
tion, Icn{r, E), is simply the integral of qcn(r, E) along the 
line-of-sight. This prescription is derived from the average 
emission profile for a sample of simulated clusters for a real- 
istic choice of parameter values (e.g., for the maximum shock 
acceleration efficiency, C, p , m ax)- In addition to the uncertain- 
ties in the model parameters there is also uncertainty in the 
observationally derived halo mass and gas density profile. 
In this work, we simply assume that the shape of qcR{r, E) 
is given by the model described above and account for the 
uncertainty in the model parameters, as well as sample vari- 
ance with an additional normalization parameter, cxcr, so 
that, 



Icr{t,E) = a C R 



qcn(r,E) 

4% 



dl. 



(11) 



We take acR = 1 as our fiducial CR model and also con- 
sider the case when ocr is fitted from the actual gamma-ray 
data as an optimal model. In the right panel of Fig. [3] we 
compare the CR profile for the fiducial model to the ex- 
pected DM annihilation profile within our three clusters, 
assuming a fiducial DM particle model with particle mass, 
M w lOOGeV, annihilating through the bb channel with 
cross-section, < av >= 3 x 10 _26 cm 3 s _1 . In general the 
CR emission is more centrally concentrated than the anni- 
hilation profile since the CR trace the gas profile. It can 
be readily seen that Fornax has a particularly low CR level 
while Coma is CR dominated. Coma has steeper profiles due 
to its larger distance and hence smaller angular size. 



10 http://www.darksusy.org. 

11 http://home.thep.lu.se/ torbjorn/Pythia.html 
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Figure 2. Photon yields for bb (left), /x+ fi~ (middle) and t+t — (right) channels. Plotted are the total photon yields including continuum 
secondary photons, final state radiation from charged final state particles, as well as inverse Compton scattering of CMB photons by 
electron/positron pairs, for the mass range 10 — 1000 GeV for the bb channel, lGeV — lOTeV for the /z+ fi~ channel and 2GeV — lOTeV 
for the t + t — channel. The masses are sampled uniformly in a log scale. Note that each spectrum cuts off at an energy corresponding 
to the particle mass. For comparison, the black line in each panel shows the photon spectrum from cosmic ray induced photons with 
arbitrary normalization. 




Figure 3. Cluster photon profiles. Left: theoretical and PSF-convolved J profile for Virgo. The total annihilation profile is shown as a 
black solid line and is decomposed into the smooth main halo part (red dashed line) and the subhalo part (blue dashed line). The green 
solid line shows the annihilation profile after PSF convolution, plotted down to an inner radius comparable to the pixel size of 0.1 deg. 
Right: PSF-convolved photon profiles from annihilation (solid) and cosmic rays (dashed) for three clusters (indicated by different colours). 
Solid lines show the expected photon count profile for the extended DM annihilation model. Dashed lines show the expected cosmic-ray 
induced photon counts for the fiducial CR model. For comparison, we also plot the PSF-convolved profile for a central point source 
(black solid line) with arbitrary normalization. In both panels, a dark matter model with particle mass, M fti lOOGeV, and annihilation 
cross-section, < av >= 3 X 10 _26 cm 3 s _1 , through the bb channel is assumed. The PSF convolutions are done with the gtmodel tool in 
the Fermi Science Tools software package. 
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3 RESULTS 

3.1 Constraints on CR emission 

With all the model components defined above, we first pro- 
ceed with ML fitting for a model with no DM annihilation 
but with cosmic rays, the "CR-only" model hereafter. Note 
that the GAL and EG backgrounds, as well as the nearby 
point sources, are always included in the analysis, as de- 
scribed in section [l~2"l The results for the CR-only model fits 
are listed in Table [2] The fitted CR levels all agree within a 
factor of three with the theoretical predictions. While For- 
nax is most consistent with no CR emission due to its in- 
trinsically low CR level, the derived upper limit for Coma 
already rules out the fiducial value at 95% confidence. 



3.2 Constraints on DM annihilation 

Given the low significance of the CR detection in the CR- 
only model, it is not safe simply to adopt the best fit olcr 
values for further extraction of the DM signal. Instead, we 
consider the following four families of cosmic ray models in 
the presence of a DM component: 

Fiducial-CR model. The CR level is fixed to the theoreti- 
cal expectation, ctcn = 1. Since this value exceeds our de- 
rived upper limit for Coma, we exclude Coma from further 
discussion of this family. 

Optimal-CR model. The CR level is taken as the best-fit 
value listed in Tabled 

Free-CR model. The normalization of the CR level is left 
as a free parameter in the fit. 

No-CR model. No CR emission is considered, only DM. 

For each family, both point source (PT) and extended 
(EXT) profiles are considered for the DM component (the 
former merely for comparison with earlier work). Note that 
when calculating the TS for DM, the null hypothesis refers 
to the full model excluding only the DM component, or 
equivalently, to the base model plus a CR component mod- 
elled according to one of our four families of CR models. We 
show results for the bb, fi + n~ and t + t~ DM annihilation 
channels. 

For none of the combinations of DM and CR models 
considered here, do we obtain a detection of DM at high 
significance in any of the three clusters. The highest sig- 
nificance is obtained for Virgo for the bb channel in a DM 
model that has a particle mass of 30 GeV and the EXT pro- 
file, in the absence of CR. In this case, we find TS = 11.6, 
corresponding to 3.4<r. This reduces to 2.6a in the Free-CR 
model and to less than la in the Fiducial-CR model. 

The value of TS =11.6 for the no-CR model for Virgo 
can be compared with the value of TS = 24 reported in 
an earlier version of this paper (arXiv:1201.1003v.l) from a 
similar analysis of the 2-year Virgo data (see Fig.0|. The de- 
crease in significance is entirely due to the subtraction of the 
new point sources which we have detected in the Virgo re- 
gion and which were not catalogued in the 2FGL. These pre- 
viously undetected sources happen to lie within the virial ra- 
dius of Virgo and can mimic the extended emission expected 
from DM annihilation. In fact, fits assuming an EXT profile 
but a power-law, rather than a DM annihilation spectrum, 
result in a similarly high significance detection, TS = 21, 
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Figure 4. The significance of a DM component in Virgo, with 
bb final states, in the absence of CR. The solid line shows the 
TS when only 2FGL sources are included in the model, while the 
dashed line shows the case when the four new point sources that 
we have detected are also included. 

and a best-fit spectral index F = —1.9. This is the typical 
spectral index of Fermi point sources (including the newly 
detected ones). The preference for a 30 GeV DM particle 
mass in the DM fits reflects a preference for a T — —1.9 
spectrum around 1 GeV, the energy scale from which most 
of the significance arises. 

In fact, the significance of the Virgo detection is further 
reduced when we take account in the analysis of a possi- 
ble undetected point source population, as we shall do in 
Appendix [B] Thus, in what follows we use our analysis ex- 
clusively to set upper limits on the flux and annihilation 
cross-section. 

8.2.1 The bb channel 

In Fig. [S] we show the 95% confidence upper limits on the 
DM annihilation flux and compare them to the CR levels. 
For each cluster, the coloured stripes are defined by the min- 
imum and maximum upper limits corresponding to the four 
families of CR models. The optimal CR levels in the three 
clusters are all comparable to the fitted DM flux, and the 
DM flux upper limits for the four different CR models vary 
only within a factor of two, with the No-CR and Fiducial- 
CR cases predicting the highest and lowest upper limits. 
The left and right panels show the results for the EXT and 
PT models respectively; the PT models always have lower 
flux upper limits than the extended models. 

The flux upper limits are translated into cross-section 
upper limits in Fig. |S] using Eqn. [3] These are also shown 
as coloured regions reflecting the variation in the different 

12 In Coma, where the Fiducial-CR model is ruled out, the 
Optimal-CR model yields the lowest upper limit. 
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Table 2. Fits to the CR-only Model 







a-CR,UL h 


F C R,UL%ph ■ cm 2 s 1 ) 


TS 


TO d 

J - ^corrected 


a CR,UL, corrected 


Coma 


0.3 ±0.1 


0.5 


2.4e-09 


5.2 


2.6 


0.6 


Fornax 


0.9 ±2 


4.8 


1.8e-09 


0.2 


0.1 


6.4 


Virgo 


0.6 ±0.3 


1.2 


2.1c-08 


8.4 


2.8 


1.6 



"Best fit normalization {acRJit = 1 is the theoretical prediction) 
'95% upper limit (UL) on the normalization 
c 95% upper limit on the CR induced gamma-ray flux from 100 MeV to 100 GeV 
d TS after allowing for undetected point sources; see Section 13.31 for details 
e Upper limit on the normalization factor after allowing for undetected point sources; see Section [33] for details 




M x /GeV M x /GeV 

Figure 5. DM annihilation flux upper limits for the bb channel. The stripes are defined by the minimum and maximum upper limits 
given by the four CR model families, with different colours corresponding to different clusters, as indicated in the legend. Left and right 
panels are the results for the EXT and PT profiles respectively. For each cluster, a solid line of the corresponding colour shows the 
optimal CR flux. 



treatments of CR. Although the predicted flux upper limits 
decrease slowly with DM particle mass and remain within 
the same order of magnitude for the mass range considered, 
the resulting cross-section upper limits increase by a factor 
of 100 from low to high particle mass. This is because low 
mass particles correspond to higher DM number densities 
(the p^/M 2 factor in Eqn.|3} for a given mass density, so to 
obtain the same flux level, the required cross-section must be 
smaller for low mass particles. With an enhancement of or- 
der 10 3 due to subhalos, a much lower cross-section is needed 
(by a factor of at least 100) for extended annihilation mod- 
els to achieve a slightly higher flux upper limit than point 
source models. 

Our cross-section limits drop below the fiducial thermal 
cross-section of 3 x 10" 26 cm 3 s _1 for M x ^100GeV. Of the 
three clusters, Virgo has the highest flux upper limits but 
it still places the tightest constraints on the annihilation 



cross-section. Our limits are muc h lower than those in th e 
11-month Fermi-LAT analysis by [Ackcrma nn et alj (2010), 
where the tightest constraint came from Fornax for a much 
lower assumed subhalo contribution of ~ 10. Our limits are 
also tighter than that from a joint analysis of the dwarf satel- 
lites o f the Milky Way by iGeringer-Sameth fc Koushiappasl 
Note that the difference between our results and 
that o f lAckermann et al.l (|201Ch comes mostly from different 
assumptions about the effect of subhalos, and only secondar- 
ily from the larger amoun t of data we have analysed. Also, 
note that in the analysis of lGeringer-Sameth fc Koushiappasl 



1 If syste matic uncertainties in the halo mass pa rameters as- 
sumed by Gcringcr-Samcth & Koushiappa^ d201 if) are consid- 
ered, the lower bounds of their derived limits become comparable 
to our limits. 
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M x /GeV 

Figure 6. Upper limits for the DM annihilation cross-section in the 
bb channel. The different colours represent the three clusters, with 
the stripes spanning the range between the minimum and maximum 
upper limits given by the four different ways of treating the CR 
component. The three highest stripes show the PT model constraints 
and the three lowest the EXT model constraints. We also plot with 
dashed lines c onstraint^! from a joint analysis of th e Milky Way 
dwarf galaxies ( Gcringcr-Samcth & Koushiappas 2011, black dashed 
line) an d previous constraints fro m the 11-month Fermi-LAT data for 
Fornax (Ackcrmann ct al. 2010, green dashed line) assuming these 
authors' optimistic value for the total enhancement due to subhalos, 
which gives the tightest constraint. The black solid line indicates the 
canonical thermal cross-section of 3 X 10 — 26 cm 3 s _1 . 



a The "Fermi- lyr" constraint is only reproduced schematically, by 
reading out several data points from the original plot in the reference. 
^ The "Fermi- lyr" constraint is only reproduced schematically, by 
reading out several data points from the original plot in the reference. 

(2011), no boost from subhalos within the halo of dwarf 
galaxies was assumed. 



3.2.2 The channel 

As have been seen in section 13.2.11 the EXT model places 
tighter constraints on the cross-section than the PT model, 
and is the fiducial model expected from recent simulations. 
Therefore from now on we will only show results for the EXT 
model. The flux and cross-section upper limits for DM an- 
nihilating through the /i + fjT channel are plotted in Fig. s[7] 
and [8] The predicted flux upper limits for Coma and Virgo 
are still comparable to the CR level, with Fornax having 
much lower CR emission. The inferred cross-section falls be- 
low the canonical value for DM particle masses less than 
10 GeV. Note the discontinuity in the upper limits around 
100 GeV which reflect the transition from the prompt an- 
nihilation dominated regime to the IC emission dominated 
regime in the photon spectrum. 



10" 7 l M 




M x /GeV 

Figure 7. DM annihilation flux upper limits in the fi + fjT channel 
for the EXT model. Line styles are as in Fig. [5] 



1 \ . > H 




M x /GeV 

Figure 8. Upper limits for the DM annihilation cross-section 
in the fi + fj,~ channel. Line styles are as in Fig. H but 
only the EXT results are s hown. The green dash ed line is 
the 11-month Fermi result jAckermann et all l2010h for For- 
nax wh ile the black dashed line is the dw arf galaxy con- 
straint (Gcrinecr-Samcth & Koushiappas l201lh . both for the 
fi + fj.~ channel. 



Extended Gamma-ray Emission in Clusters 11 



10 



10 



10 



-23 



-24 



-25 



s 10 

b : 

V 



10 



10 



-27 



-281 



I Coma 
Fornax 
I Virgo 
■ Dwarf 



10' 



10' 

M x /GeV 



10° 



10 H 



Figure 9. Upper limits for the DM annihilation cross-section 
in the t+t — channel. Line styles are as in Fig. \E\ but only the 
EXT result s are shown. The black dashed line is the dwarf galaxy 
constraint. l|Geringer-Sameth fc Koushiappa s 2011) 



3.2.3 The t+t~ channel 

In Fig. [9] we show the cross-section upper limits for the 
t + t~ channel. This i s the primary component of the lep- 
tonic model used by iHooper fc Linden] l|2011bT ) to fit the 
excess gamma-ray emission in the Galactic center region. 



3.3 Allowing for an undetected point source 
population 

Although we have detected five new point sources in the 
45-month data in the region of our three clusters, it is still 
necessary to account for population of still undetected point 
sources. When no unknown point sources are present, the 
probability of measuring a certain value of TS when the 
null hypothesis is true is given by the probability that Pois- 
son fluctuations in the photon counts for the null model 
exceed some value. When a population of undetected point 
sources is present, the Poisson fluctuations become corre- 
lated and it is easier for the same amplitude of fluctuations 
to result in a given value of TS. In this case, the distribution 
of TS no longer follows the x 2 distribution, as predicted by 
Wilk's theorem, because the data no longer follow a pure 
Poisson distribution from which the likelihood function is 
constructed. 

Allowing for the presence of undetected point sources in 
the data will lead to weaker upper limits. We obtain these by 
performing Monte-Carlo simulations to re-calibrate the sig- 
nificance corresponding to a given value of TS. In the simu- 
lations we include the GAL and 2FGL sources, but we split 
the EG component into two parts: a population of unde- 
tected point sources and a residual smooth EG component, 
such that the sum of the two is consistent with the standard 



EG component. We consider a benchmark model for the un- 
detected p oint source populat ion which is close to the model 
derived bv lAbdo et al.l l|2010t) and which contributes 14% of 
the EG background. Details of the simulations may be found 
in Appendix [B] A standard likelihood analysis is then per- 
formed on the simulated data in order to derive appropriate 
values of TS for an assumed DM or CR component. 

With the introduction of the undetected point source 
population, the distribution function for TS is found to be 
roughly described by '{TS /b) /2. That is, the significance 
of a given value of TS is approximately reduced by a factor 
of b compared to the significance of the same value of TS 
in the absence of the undetected point source population. 
For a DM component, we find b ~ 2 for Coma and b ~ 3 
for Virgo and Fornax. The b factor is not sensentive to the 
adopted DM spectrum. For the CR models, we find b ~ 2 
for Coma and Fornax, and b ~ 3 for Virgo. 

In order to obtain new limits from the corrected TS, let 
us first consider the likelihood function that has been max- 
imized over all the nuisance parameters. Expanding around 
the maximum likelihood value of the parameter N to leading 
order, we have 



lnL(iVo) = hiL(A0 - -H{N 



N) 



(12) 



where H = 



df In L 

dJV 2 



is the Hessian maxtrix. The 95% 
InL(N) = -1.35, 



upper limit is calculated from \nL(UL) 
so that 

UL = N + 1.64a N . (13) 

Note that — = VTS. Similar equations hold for the inl- 
aw 

proved background model, with UL, on and TS replaced 
by UL', a' N and TS' respectively, assuming that there is 
no bias in the best-fit parameters. The 95% upper limit is 
then corrected for the undetected point source fluctuations 
according to 



UL' = 



/TS+lMy/TS'/TS 



UL, 



(14) 



VTS+1M 

where UL and TS are the upper limit and likelihood ra- 
tio from the standard analysis, while UL' and TS' are the 
corrected upper limit and likelihood ratio. For b = 3, the 
increase in the upper limit is at most 70%. 

In Fig. 1101 we show the corrected dark matter annihila- 
tion cross-section upper limits adopting b = 2 for Coma and 
b = 3 for Virgo and Fornax. The corrected TS and upper 
limits for CR models are listed in the last two columns of 
Table [2] 
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Figure 10. Upper limits for the DM annihilation cross-section in the bb (left), (middle), and t + t~ (right) channels, after 

including the effect of undetected point sources. Line styles are as in Fig. [6] but only the EXT results are shown. Note that the lower 
bounds of each band are still determined by the results without including undetected point sources in the analysis. 



4 DISCUSSION AND CONCLUSIONS 

We have performed maximum likelihood fits to the 3-year 
Fermi-LAT data for three galaxy clusters: Coma, Fornax and 
Virgo. We fit models which, in addition to point sources 
and galactic and extragalactic backgrounds, include emis- 
sion due to dark matter (DM) annihilation and cosmic rays 
(CR). For the former, we assume both a point source and the 
theoretically predicted extended distribution of gamma rays 
in three generic annihilation channels, the bb, fj> + fj,~ and 
t + t~ channels. When searching for a dark matter signal, 
we experiment with different treatments of the CR com- 
ponent. In the traditional Fermi analysis, the extragalactic 
background (EG) is assumed to be a smooth component. 
In this work we have also investigated a more realistic EG 
model where a fraction of the EG emission comes from a 
population of undetected point sources. 

Performing a standard likelihood analysis we obtain the 
following results: 



Way dwarf galaxies llGeringer-Sameth fc Koushiappasll201ll ; 
iThe Fermi-LAT Collaboration: M. Ackermann et al.ll201ll ). 

Our new limits exclude the thermal cross-section for M x < 
100 GeV for bb and t + t~ final states, and for M x < 10 GeV 
for fi~ final states. We note that the annihilation cross- 
section in dark matter halos need not be the standard ther- 
mal cross-section of supersymetric models. In cases where 
the cross-section is velocity dependent, for e xample, through 
p-wav e contributions at freeze-out (see e.g.. lJungman et alj 
( 1996)), one can easily have a different average cross-section. 
The thermal cross-section, however, could still be reconciled 
with the data by assuming a larger cutoff mass in the WIMP 
power spectrum, thus reducing the contribution from sub- 
halos and hence the J factor. Since the tota l enhancement 
from subhalo emission scales as b oc M~^ t 226 (|Springel et alj 
120081 ), a cut-off mass of 1O _4 M0, rather than our assumed 
10 _6 Mq, would be sufficient to increase the cross-section 
limits by a factor of 3. 



(i) In all three clusters and for the four different treat- 
ments of CR we have implemented, no significant detection 
of DM emission is obtained. We set upper limits on the flux 
and cross-section of DM annihilation in the three clusters 
we have investigated. Uncertainties in the CR component 
have only a mild effect on the upper limits: for the different 
CR models, the DM upper limit constraints agree to within 
a factor of two. 

Models in which the DM annihilation emission has 
the e xtended profile predicted by cosmological simula- 
tions (|Gao et al.ll2Qlll ) have higher flux upper limits than 
models in which this emission is assumed to be a point 
source. Due to the large luminosity enhancement, of or- 
der of 1000, by emission from subhalos, the upper lim- 
its on the annihilation cross-section for extended models 
are at least 100 times lower than those for point source 
models. Our cross-section constraints are much tighter 
than those f rom an analysis of cl usters using the 11- 
month data (|Ackermann et all 12010] ). mostly because we 
take into account the effect of subhalos. Our constraints 
are also tighter than those from a joint analysis of Milky 



(ii) Assuming no DM annihilation radiation, the gamma 
ray data for Coma and Virgo already set significant con- 
straints on the CR level. For Virgo, the data are consistent 
with the predictions of th e analyt ic CR mode l prop osed by 
iPinzke fc Pfrommeij (|201Ch and IPinzke et all (|201ll ) while, 
for Coma, the data place an upper limit that is a factor of 
two below the analytical prediction, indicating either an un- 
certainty in model parameters such as halo mass, gas density 
and maximum shock injection efficiency, £ p , m ax, or a pecu- 
liarity of the CR emission in Coma. If attributed to £ P: ma,m> 
the upper limit on the normalization parameter, acn, trans- 
lates into an upper limit on Q Pt m,ax of 0.3, assuming a lin- 
ear form for g(C, Pt max). This is consistent with th e estimates 
obtained independently by IZimmer et a l. (2011) for Coma 
using Fermi data and by the IMAGIC Collaboration et al.l 
(|20 1 IT ) for the Persus cluster using MAGIC observations. If 
interpreted as an error in the halo mass, a reduction in mass 
by a factor of 1.6 is required to reconcile the model with the 
upper limits, assumi ng a simple CR luminosit y scaling re- 
lation, Lj oc M200 6 (|Pinzke fc Pfrommerll2010l ). or a factor 



of 4.3 according to Eqn. ID1I in the case when the gas den- 
sity profile is fixed from X-ray observations. For Fornax, the 
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zero-significance of a CR component is consistent with the 
low level predicted by the model. 

(iii) Five new point sources with TS > 25 in Virgo and 
Fornax have been detected in the 45-month data. Ignoring 
these new point sources results in a ~ 5a detection for a DM 
component in Virgo, in contrast to a ~ 3<r detection when 
account is taken of these point sources. 

In addition to the standard likelihood analysis, we 
have also investigated a model in which the EG component 
includes a population of undetected point sources whose 
number-flux relation extrapolates smoothly that of the de- 
tected sources. Using Monte-Carlo simulations, we find that 
the standard Fermi likelihood analysis could overestimate 
the TS of extended emission by a factor of 2 — 3, and un- 
derestimate the upper limits by up to 70 percent. Adopting 
this more realistic EG model yields slightly looser upper 
limits, but does not quantitatively change any of the above 
conclusions. Still, it should be kept in mind that these cor- 
rections are derived from simulations assuming a particular 
distribution for the point source population. It is too compu- 
tationally expensive to explore the parameter space of point 
source populations with Monte-Carlo simulations. A more 
detailed and more general analytical study of the effe ct of 
undetect ed point sources will be presented elsewhere |Hanl 
I in prep] ). 

In our analysis we have allowed the parameters of 2FGL 
point sources lying within the cluster virial radius to vary. 
This accounts for possible corrections to the 2FGL parame- 
ters in the presence of a DM or a CR component, while also 
avoiding the risk of refitting sources lying near the bound- 
ary of the data region with less accuracy. The parameters 
of highly variable sources are also kept free since the 2FGL 
parameters for these sources would be the average during 
a 2 year period whereas here we have 45 months of data. 
However, we also tried keeping all the point sources fixed or 
allowing the parameters of all the point sources within the 
data region to vary during the fitting. We find that this free- 
dom in the treatment of the point sources has little impact 
on the DM model fits. 

The cluster annihilation luminosity scales roughly lin- 
early with halo mass, with the shape of the profile being 
almost independent of halo mass or concentration when ex- 
pressed in terms of the normalized radius r/i?200- We inves- 
tigate the effect of mass uncertainties in Appendix [E] We 
have also checked that th e different energy cu ts assumed in 
our analysis and in that of lHuang et al.l (|201 lh have no effect 
on the derived upper limits. We are able to repr oduce the up- 
per lim its on the annihilation cross-section of lHuang et al.l 
|201l[ ) for the test case of the Fornax cluster with 3-year 
data, after adopting the same instrument response function 
and correcting for slightly different assumed subhalo contri- 
butions. 

The CR model used in this analysis is still subject to im- 
provement. This model is derived from simulations which, 
unavoidably, make simplifying assumptions. For example, 
the simulations only consider advective transport of CR by 
turbulent gas motions but there are other processes such as 
CR diffusion and stream ing which may flatten the CR pro- 
files (|Enfilin et al.ll201lh . In particular, if the CR diffusion 
is momentum dependent this will entangle the spectral and 
spatial profile of CR and modify the morphology as well as 



the spectrum of the CR emission, thus invalidating our basic 
assumption that acR is the only free parameter. There could 
also be CR injected from AGN which are not accounted for 
in the current model. 

Although we have not detected DM annihilation emis- 
sion in our small cluster sample, the signal-to-noise ratio can 
potentially be enhanced by stackin g many clusters. Su ch an 
analysis was recently carried out bv lHuang et al.l (|201ll ). but 
the signal-to-noise was degraded because of their assumption 
of an NFW annihilation profile. These authors considered an 
extended subhalo-dominated annihilation profile but only 
for individual clusters, not for the stack. Their stacked anal- 
ysis placed looser constraints on DM annihilation emission 
than their analysis of individual clusters, presumably be- 
cause the use of an inappropriate theoretical profile resulted 
in the different clusters yielding inconsistent results. Thus, 
it is clearly worth repeating the joint analysis with the "cor- 
rect" subhalo-dominated profile. It is also tempting to ex- 
tend the search for DM annihilation using multi-wavelength 
data, from the radio to very high energy gamma-rays, where 
different systematics are expected for different bands. 
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APPENDIX A: DETECTION OF NEW POINT 
SOURCES 

We model the new point sources assuming power-law spec- 
tra. For a given pixel, we calculate the TS value for an as- 
sumed new point source centered on that pixel. The TS 
calculation is performed using the binned method in the 
pyLikelihood tool, with a null model which includes the 
GAL and EG components and all the 2FGL sources within 
15 deg of each cluster, but with the parameters of the 2FGL 
sources fixed. Around each cluster, we carry out a first scan 
of all the pixels within the cluster virial radius (and within 
4 deg around Coma) using a pixel size of 0.2 deg. 

Regions with a peak TS > 16 are identified as potential 
locations of new point sources. We then scan each poten- 
tial point source region using 10 times smaller pixels. The 
calculated TS map is then interpolated with cubic splines 
down to 0.002 deg/pixel. The value and location of the TS 
peak is taken as the TS and position for a new point source, 
if the peak TS > 25. In case several peaks are clustered, 
we first extract the primary TS peak, then scan for lower 
TS peaks by including the newly detected sources into the 
null model. In our sample, no secondary peaks survive this 
iterative examination to be identified as new point sources. 

The new point sources are listed in Table lATI and ploted 
in FigurefTland lCll Sources in Virgo and Fornax are prefixed 
by "V" and "F" in their names respectively. None of these 
new sources show significant variability when binned over 
monthly scale. The last column of Table IA1I shows possible 
associations of astrophysical sources with these new detec- 
tions, which are found to lie within the 2a confidence region 
of the detections. 



Extended Gamma-ray Emission in Clusters 15 



Table Al. Newly detected point sources 



Name 


TS 


RA (deg) 


DEC (deg) 


Flux (10~ 9 ph ■ cm" 2 s" 1 ) 


Spectral Index a 


Seperation (deg) b 


Possible Association 


VI 


32.5 


190.920 


16.194 


5.9 ± 1.4 


-2.3 ±0.2 


4.96 


LBQS 1241+1624 


V2 


31.8 


185.698 


11.116 


3.7 ± 1.0 


-2.0 ±0.2 


2.31 


[VV2006] J122307.2+110038 


V3 


31.6 


184.066 


9.456 


2.3 ±0.8 


-1.9 ±0.2 


4.58 


2MASX J12160619±0929096 


V4 


30.5 


185.894 


8.286 


1.6 ± 0.7 


-1.8 ±0.2 


4.42 


SDSS J122321.38+081435.2 


Fl 


26.3 


58.300 


-36.386 


0.9 ±0.6 


-1.7 ±0.3 


3.17 


[VV98b] J035305. 1-362308 



a Photon spectral index /3 for dN/dE oc . 
''Distance to cluster center. 



APPENDIX B: MONTE-CARLO SIMULATION 
OF UNDETECTED POINT SOURCE 
POPULATIONS 

To model the undetected po pulation we adop t the following 
model based on the results of lAbdo et all {2010). Each point 
source is assumed to have a power-law spectrum defined by 
two parameters: flux and spectral index. The spectral index 
distribution is modeled as a Gaussian of mean /i = 2.36 and 
a = 0.27. The flux and spectral index are assumed to be 
independent. The differential number density of undetected 
point sources is assumed to be given by 

We adopt S b = 6.6 x 10" 8 ph cm -2 s"\ /3 = 1.58 and 
A = 4.1 x 10 8 c m 2 s Sr -1 , as derived from Table 4 of 
lAbdo et al.l (|2010t ). Since the total number of point sources 
diverges for ft > 1, we cut off the flux distribution at 
Smin = 1 x 1CP 11 ph cm -2 s _1 . Due to the dependence of 
the detection efficiency on flux and spectral shape, there 
is no obvious cutoff in the maximum flux of undetected 
sources. We take S m3x = 1 x 10 -8 ph cm -2 s _1 as the 
detection threshold which corresponds to a detection com- 
pleteness of ~ 50%, comparing 2FGL source counts and the 
model. This implies an undetected point source flux of 14% 
of the standard EG b ackground, consistent with the results 
of lAbdo et ail l|2010h . The synthetic spectrum of these un- 
detected point sources is then subtracted from the standard 
EG template to yield a residual EG template for the simu- 
lation. 

We perform 750 independent realizations of the 15 deg 
Virgo region in the presence of undected point sources. For 
each realization, we generate mock data in the following 
steps: 

• Generate a Poisson random number for the total num- 
ber of undetected point sources within 15 deg. 

• For each point source, generate a random spectral in- 
dex and a random flux according to the distributions speci- 
fied above. Also, generate random coordinates for the point 
source according to a uniform distribution on the sky. 

• Feed these point sources and the 2FGL point sources 
within 15 deg, as well as the GAL and remaining EG com- 
ponents, to gtobssim. 



simulated data without including any of the randomly gen- 
erated point sources in the model. Here we only consider 
the CR-only and the DM-only models. In Fig. IB1I we show 
the cumulative probability distribution of TS values. Sim- 
ple scaled versions of the standard x 2 (TS) /2 distributions 
can roughly describe the TS distribution and provide the 
simplest way to convert the fitted TS to the standard \ 2 - 
distributed TS. 



The standard likelihood analysis is then applied to the 
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Figure Bl. Distribution of TS from simulated data which include a population of undetected point sources. Left: the distribution of 
TS for DM-only models, where the DM particle mass is taken to be ~ 30 GeV and the DM follows the EXT cluster profile. Right: the 
distribution of TS for CR-only models. In each panel the dashed lines show the distribution extracted from the simulations for three 
cluster models and the solid lines show a rescaled version of the standard cumulative \ 2 distribution. 



APPENDIX C: GAMMA-RAY IMAGES FOR 
COMA AND FORNAX 

In this Appendix we show gamma-ray images for the Coma 
and Fornax cluster regions. The corresponding image for 
Virgo is shown in Fig. [1] 
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Figure CI. Integrated gamma-ray images in the Coma (left) and Fornax (right) cluster regions. The green dashed circle marks the 
virial radius of the cluster. Each image covers 20 X 20 deg 2 with a pixel size of 0.1 deg, constructed from the 3-year Fermi-LAT data 
applying the data cuts described in the main text. The small solid circle in Fornax marks the position of a newly detected point source. 



APPENDIX D: SEMI- ANALYTIC FORMULA 
FOR THE COSMIC RAY INDUCED 
GAMMA-RAY EMISSION 

Here we summarize the relevant equations for calculating 
the CR indu ced gamma-ray emi s sion in ga laxy clusters 
as der ived bv lPinzke fc Pfrommerl (|20T0t ) and IPinzke et all 
|201ll ). The CR induced photon source function from pion 
decay can be decomposed as: 



A(r)s(E). 



dtdVdE 
The spatial part is given by: 

A(r) = ((C 2 oo-C' cen i er )(H-(— — — )~ fl )~ 1 +C center )pgas(r) 2 , 

-t^trans 

(Dl) 

with 

Ccenter = 5 X 10" 7 (D2) 

C200 = 1-7 x 1(T 7 x (M 20 o/10 15 M Q ) - 51 (D3) 
Rtrans = 0.021i? 20 o x (M 200 /10 15 M ) °' 39 (D4) 
P = 1.04 x (M 2O o/lO 15 M ) 015 (D5) 

The spectrum is given as: 

= g((^p,max)D^(E^ , Ey : break) 3m ij ~ 

^3 o- 



tral pion mass and c the speed of light. The maximum shock 
acceleration efficiency is chosen to be ( p , ma x = 0.5 so that 
g(C P ,max) = 1. The term D 7 (_E 7 , E^^ reak ) describes the dif- 
fusive CR losses due to escaping protons as 



Dj(Ej, Ejfireak) — [1 + ( 



The proton cut-off energy is 



E 



-,31-1/9 



E- 



7, break 



E pM ~ _GeV( T ^) . 



(D7) 



(D8) 



The energy Ep^r&ak is related to the photon cut-off energy, 
E-f, break, through the momentum relation P 7 m The 
effective cross-section for proton-proton interactions is given 
by: 

a PP}l ~ 32(0.96 + e 4A2 - 2A °" ) m barn. (D9) 
The gas density is fitted with multiple beta-profiles: 



Pgas 



? { J2nU0)[l + ( — ) 2 ]- 3P <} 1/2 , (D10) 



with A = (0.767,0.143,0.0975), a = (2.55,2.3,2.15), <5< ~ 
0.14cv~ 1.6 + 0.44. Here m p is the proton mass, m n o the neu- 



where Xh = 0.76 is the primordial hydrogen mass fraction 
and X e = 1.157 is the ratio of electron and hydrogen num- 
ber densities in the fully ionized intracluster medium, with 
p arameter values for Wj(0), r c ,i and /3j listed in TABLE VI 
of IPinzke et"afl l|201ll ). 
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Figure El. Upper limits on the cross-section for DM annihilating 
into the bb channel in the no-CR model. Different lines correspond 
to different adopted values for the mass of the dark matter halo of 
the Virgo cluster, as labelled. No allowance for undetected point 
sources has been made in this figure. 



APPENDIX E: EFFECT OF MASS 
UNCERTAINTIES IN VIRGO 

We adopt a viria l mass for Virg o of 7. 5 ± 1.5 x 1O 14 M0, 
as estimated by iTullv fc Shavaf dl984h from an analysis 
of the infall pattern of galaxies around Virgo. Thi s value 
is consistent with other dynamical measurements (jSmithl 
19361 ; iHoffman et al.lll980l; [Tonrv et al.ll200Cl : iFouque et alj 



200ll ; lKarachentsev fc Nasonovall201fj| '). Mass estimates from 



X-ray gas modelling t e nd to give somewhat lower values 
jBohringer et all Il994l ; ISchindler etHI Il999l ; | Urban et al.l 
l201ll ). Thus, in add i tion t o our adopted mass uncertainty 
from Tullv fc Shaval(|l984h . as an extreme case, we consider 
also a value of 1.4 x 1O 14 M0, obtained by scalin g the X-ray 
estimate to the virial radius (|Urban et al.ll201lh . In Fig. IE1I 
we show the effect of adopting these different masses on the 
upper limits for DM annihilation in the bb channel. Since the 
flux upper limit is insensitive to slight changes in the profile 
shape and thus in the mass, while the luminosity (or inte- 
grated J factor) scales linearly with mass, the cross-section 
upper limits are expected to be roughly inversely propo- 
tional to mass. This is indeed the case in Fig. IE1I where 
< av >ul<x M^° 9 . 



